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Efficient polynomial time algorithms computing 
industrial-strength primitive roots 

Jacques Dubrois*and Jean-Guillaume Dumas ^ 
December 9, 2008 



E. Bach, following an idea of T. Itoh, has shown how to build a small 
set of numbers modulo a prime p such that at least one element of this set 
is a generator of Z/pZ. E. Bach suggests also that at least half of his set 
should be generators. We show here that a slight variant of this set can in- 
deed be made to contain a ratio of primitive roots as close to 1 as necessary. 



algorithm providing primitive roots of p with probability of correctness 
greater than 1 — e and several 0{log" [p)), a < 5.23, algorithms comput- 
ing "Industrial-strength" primitive roots. 

1 Introduction 

Primitive roots are generators of the multiplicative group of the invertibles of a 
finite field. We focus in this paper only on prime finite fields, but the proposed 
algorithms can work over extension fields or other multiplicative groups. 

Primitive roots are of intrinsic use e.g. for secret key exchange (Diffie- 
Hellman), pseudo random generators (Blum-Micali) or primality certification. 
The classical method of generation of such generators is by trial, test and 
error. Indeed within a prime field with p elements they are quite numerous 
{4){4>{p)) = 4>{p — 1) among p — 1 invertibles are generators. 

The problem resides in the test to decide whether a number g is a generator 
or not. The first idea is to test every for i = l..p — 1 looking for matches. 
Unfortunately this is exponential in the size of p. An acceleration is then to 
factor p — 1 and test whether one of the g i is 1 for q a divisor of p — 1. If this 
is the case then g is obviously not a generator. On the contrary, one has proved 
that the only possible order of is p — 1. Unfortunately again, factorization is 
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still not a polynomial time process: no polynomial time algorithm computing 
primitive roots is known. 

However, there exists polynomial time methods isolating a polynomial size 
set of numbers containing at least one primitive root. Shoup's [23] algorithm 
is such a method. Elliot and Murata [S] also gave polynomial lower bounds on 
the least primitive root modulo p. One can also generate elements with expo- 
nentially large order even though not being primitive roots [13j . Our method is 
in between those two approaches. 

As reported by Bach [5], Itoh's breakthrough was to use only a partial fac- 
torization of p — 1 to produce primitive roots with high probability |15| . Bach 
then used this idea of partial factorization to give the actually smallest known 
set, deterministically containing one primitive root [2], if the Extended Riemann 
Hypothesis is true. Moreover, he suggested that his set contained at least half 
primitive roots. 

In this paper, we propose to use a combination of Itoh's and Bach's algo- 
rithms producing a polynomial time algorithm generating primitive roots with 
a very small probability of failure (without the ERH). Such generated numbers 
will be denoted by "Industrial-strength" primitive roots. We also have a guar- 
anteed lower bound on the order of the produced elements. In this paper, we 
analyze the actual ratio of primitive roots within a variant of Bach's full set. 
As this ratio is close to 1, both in theory and even more in practice, selecting a 
random element within this set produces a fast and effective method computing 
primitive roots. 

We present in section [2] our algorithm and the main theorem counting this 
ratio. Then practical implementation details and effective ratios are discussed 
section m We conclude section |6] with applications of primitive root generation, 
accelerated by our probabilistic method. Among this applications are DifHe- 
Hellman key exchange, ElGamal cryptosystem, Blum-Micali pseudo random bit 
generation, and a new probabilistic primality test based on Lucas' determin- 
istic procedure. This test uses both the analysis of the first sections and the 
composite case. 

2 The variant of Itoh/Bach's algorithm 

The salient features of our approach when compared to Bach's are that: 

1. We partially factor, but with known lower bound on the remaining factors. 

2. We do not require the primality of the chosen elements. 

3. Random elements are drawn from the whole set of candidates instead of 
only from the first ones. 

Now, when compared to Itoh's method, we use a deterministic process producing 
a number with a very high order and which has a high probability of being 
primitive. On the contrary, Itoh selects a random element but uses a polynomial 
process to prove that this number is a primitive root with high probability [15j . 
The difference here is that we use low order terms to build higher order elements 
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Algorithm 1: Probabilistic Primitive Root 



Input: A prime p> 3 and a failure probability < e < 1. 

Output: A number, primitive root with probability greater than 1 — e. 

begin 

Compute B such that (1 + ^)(1 - -g)'"*^^ =l-e. 

Partially factor p — 1 = P'hQ (P^ ^ ^ '^'^'^ Q ^'^^ factor 

<B). 

for each 1 < i < h do 

By trial and error, randomly choose verifying: 

a.^''' ^ 1 (modp). 

h ^ 
Set a = a^^ {modp). 

i=l 

if Factorization is complete then 
I Set Probability of correctness to 1 and return a. 
else 

Refine ProbabiHty of correctness to (1 + q^)(1 — -g)'"^^ '3. 

p-i 

Randomly choose b verifying: b Q ^1 and return 

p-i 

g = ab Q (modp). 

end 



whereas Itoh discards the randomly chosen candidates and restarts all over at 
each failure. Therefore we first compute the ratio of primitive roots within the 
set. We have found afterwards that Itoh, independently and differently, proves 
quite the same within his [15, Theorem 1]. 

Theorem 1 At least of the returned values of Algorithm [7] are primitive 
roots. 

Proof. We let p — 1 = kQ. In algorithm [Tl the order of a is {p — 1)/Q — k 
(see [2]). We partition Z/pZ* by S and T where 

S ^{be Z/pZ* : h'^ ^ l{mod p)} and T = {5 e Z/pZ* : b'' = l{modp)} 

and let U — {b E Z/pZ* : b'^ has order Q}. Note that for any x G Z/pZ* of 
order n and any y G Z/pZ* of order m, if gcdin, m) = 1 then the order of 
z = xy{mod p) is nm. Thus for any b £ U it follows that g = ab^{mod p) has 
order p — 1. Since U S, we have that of the returned values of algorithm 
[T]are primitive roots. 

We thus now count the number of elements of U and S. On the one hand, 
we fix arbitrarily a primitive root g S Z/pZ* and define E = {i : < i < 
Q and gcd{i,Q) = 1}. \E\ = (p{Q) and it is not difficult to see that U = 
{g'+^Q : j e -B and < j < fc - 1}. This implies that \U\ = kip{Q). 
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On the other hand, we have T — {ff',g'^ , ■ ■ ■ ,g^^ ^^"^j. The partitioning 
therefore gives \S\ = \Z/pZ*\ - |T| = p - 1 - fc. We thus conclude that ||[ = 
fc0(Q) _ W) n 

p-l-k ~ Q-1 ■ '-' 

Corollary 2 Algorithm]^ is correct and, when Pollard's rho algorithm is used, 
has an average running time of O ^•y/^log^'^(p) + log^(p) log(log(p))^0. 



Proof. We first need to show that > 1 — e. Let Q — fl li'^^ where 

i=l 

uj{Q) is the number of distinct prime factors of Q. Then <j){Q) = Yi 4>{li^^) — 

i=l 

"(Q) '-(Q) 
Q n (1 - Thus ^ = (1 + n (1 " Now, since any factor of 

i—l i—1 

Q is bigger than B, we have: O (1 - ^) > 11 (1 - "g) = (1 - ^)'^^^'> ■ To 

i=i i=i 

conclude, we minor lu{Q) by log^^Q). This gives the probability refinemeniQ. 
Since Q is not known at the beginning, one can minor it there by since p—1 
must be even whenever p > 3. Now for the complexity. For the computation of 
B, we use a Newton-Raphson's approximation. The second step depends on the 
factorization method. Both complexities here are given by the application of 
Pollard's rho algorithm. Indeed Pollard's rho would require at worst L = 2\B~\ 
loops and L — 0{\/B) on the average thanks to the birthday paradox. Now 
each loop of Pollard's rho is a squaring and a gcd, both of complexity O(log^p). 

Then we need to bound B with respect to e. We let h = (p — l)/2 and 
B* = min{ln{h)/e; h} and consider fh{e) = (1 - l/B*y°SB*ih) _ (1 - e). Then 



hi{B*)) 2ln{B*) \ln{B*) ln{h) J \6ln{B* 

is strictly positive as soon as B* > 3. This proves that 1-e < {1-1 /B*y°SB'(h)_ 
Now, since (1 — l/_B)i°gB('') is decreasing in B, this shows that B such that 
(1 + _ ^)iogB ^ ^ 1 _ e satisfies B < B* < 

For the remaining steps, there is at worst logp distinct factors, thus logp dis- 
tinct ai, but only log log p on the average TT, Theorem 430]. Each one requires 
a modular exponentiation which can be performed with O(log^p) operations 
using recursive squaring. Now, to get a correct a^, at most O(loglogp) trials 
should be necessary (see e.g. j25l Theorem 6.18]). However, by an argument 

similar to that of theorem [U less than 1 — ^ of the are such that a^''' = 1. 
This gives an average number of trials of 1 + — , which is bounded by a constant. 



* Using fast integer arithmetic this can become : 

O logi-^{p) log2(log(p)) log(log(log(p)))+ log2(p) log2(log{p)) log{log(log(p)))) ; but the 

worst case complexity is O (i log^(p) + log*(p) log(log(p))) . 

tNote that one can dynamically refine B as more factors of p — 1 are known. 
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This gives log x log^ x log log in the worst case (distinct factors x exponentia- 
tion X number of trials) and only log log x log^ x2 on the average. □ 



3 About the number of prime divisors 

In the previous section, we have seen that the probability to get a primitive root 
out of our algorithm is greater than (l — .^j'^'''^' q ^j^g remaining unfactored 
part with no divisors less than B. The running time of the algorithm, and in 
particular its non-polynomial behavior depends on B and on to. In practice, uj is 
quite small in general. The problem is that the bound we used in the preceding 
section, logg{p — 1), is then much too large. In this section, we thus provide 
tighter probability estimates for some small B and large Q. 

Theorem 3 Let B G IN. Q ^ IN such that no prime lower than B divides Q 
then: 



Lo{Q) < \ogsiQ) VS > 2 (1) 

logB(/n(Q)) 

Proof. Of course, ([T]) is a large upper bound on the number of divisors of 
Q and therefore a bound on the number of prime divisors. Now for the other 
bounds, we refine Robin's bound on w j231 Theorem 11]: which is uj(n) < 
in{in(n)) ^"■("-)- Let Nk — Y[i=i Pi whcrc Pi is the i-th prime. Now, we let k 
be such that <Q < Then to{Q) < to [j^] =k- tt{B) since 

no prime less than B can divide Q. We then combine this with the fact that 
X '—^ lll^ is decreasing for {X > e), to get: lj{Q) < iog^^(^;^^Q)-| ^ogg{Q) where 



(fc-,r(B))log( 









log I 




1 



F{k,B) = / N ■ We then replace both Nk in F{k,B) using 



e.g. classical bounds on 9{pk) — ln{Nk) [23_, Theorems 7 & 8]: 

f, , , , , . ^" fc- 2.1454\ 

SiPk) >k { In k + ln In k - 1 -\ (5) 

\ In k J 

X , A , , , , In In fc- 1.9185\ 

e(pk)<k[lnk + lnlnk-l + — (6) 

\ In k J 

We therefore obtain a function F{k, B) explicit in k and B. The values given 
in the theorem are the numerically computed maximal values of F{k, B) as a 



5 



function of fc for B e {2^°, 2^^, 2"^^}. The claim then follows from the fact that 
F{k,B) is decreasing in B. □ 
It is noticeable that the last estimates are more interesting than log^((5) only 
when B^^'^'^^ < ln{Q). Those estimates are then only useful for very large Q 
(e.g. more than 10^ bits for B = 2^^). 



4 Industrial-strength primitive roots 



Of course, the only problem with this algorithm is that it is not polynomial. 
Indeed the partial factorization up to factors of any given size is still exponen- 
tial. This gives the non polynomial factor ^J^. Other factoring algorithms with 
better complexity could also be used, provided they can guarantee a bound on 
the unfound factors. For that reason, we propose another algorithm with an at- 
tainable number of loops for the partial factorization. Therefore, the algorithm 
is efficient and we provide experimental data showing that it also has a very 
good behavior with respect to the probabilities: 

Heuristic 2: Apply Algorithm]^ with B < log^(p) log^(log(p)). 

With Pollard's rho factoring, the algorithm has now an average bit polynomial 
complexity of : O (log^(p) log(log(p))) (just replace B by log^(p) log^(log(p)) 
and use L — ^/B). In practice, L could be chosen not higher than a million: 
in figures [1] we choose Q with known factorization and compute ; the 
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Figure 1: Actual probability of failure of Algorithm [T] with L = 2^° 
experimental data then shows that in practice no probability less than 1 — 2^^° 
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Figure 2: Actual probability of failure for Q with many distinct factors 



is possible even with L as small as 2^". 

Provided that one is ready to accept a fixed probability, further improvements 
on the asymptotic complexity can be made. Indeed, D. Knuth said "For the 
probability less than (^)^^ that such a 25-times-in-row procedures gives the wrong 
information about n. It's much more likely that our computer has dropped a bit 
in its calculations, due to hardware malfunctions or cosmic radiations, than that 
algorithm P has repeatedly guessed wrong, 'd We thus provide a version of our 
algorithm guaranteeing that the probability of incorrect answer is lower than 

2-40. 

Algorithm 3: If p is small (p < 45171967j, factor p — 1 completely, otherwise 
apply Algorithm[l\ with B = log^ '^^^^^'^p. 

With Pollard's rho factoring, the average asymptotic bit complexity is then 
OiXog^'^^'^^^'' p): Factoring numbers lower than 45171967, takes constant time. 
Now for larger primes and B — log"(p), we just remark that (1 + — 

iy°gi3 T- is increasing in p, so that it is bounded by its first value. Numerical 
approximation of a so that the latter is 1 — 2^'"' gives 5.298514. The complex- 
ity exponent follows as it is 2 + ^. One can also apply the same arguments 
e.g. for a probability 1 — 2~^^ and factoring all primes p < 2^^^ (since 513-bit 
numbers are nowadays factorizable) , then slightly degrading the complexity to 
0{log^'^^^^^^p). We have thus proved that a probability of at least 1 — 2^*° can 
always be guaranteed. In other words, our algorithm is able to efficiently pro- 
duce "industrial-strength" primitive roots. This is for instance illustrated when 

t More precisely, cosmic rays only can be responsible for 10^ software errors in 10^ chip- 
hours at sea level|20] . At IGHz, this makes 1 error every 2^^ computations. 
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Figure 3: Generations of primitive roots 

comparing our algorithm, implemented in CH — h with GMP, to existing software 
(Maple 9.5, Pari-GP, GAP 4r4 and Magma 2.1ljlon an Intel PIV 2.4GHz. This 
comparison is shown on figure [31 Of course, the comparison is not fair as other 
softwares are always factoring p — 1 completely. Still we can see the progress in 
primitive root generation that our algorithm has enabled. 



5 Analysis of the algorithm for composite num- 
bers 

In this section we propose an analysis of the behavior of the algorithm for 
composite numbers. Indeed, our algorithm can also be used to produce high, 
if not maximal, order element modulo a composite number. This analysis is 
also used section [6.21 for the probabilistic primality test. It is well known that 
there exists primitive roots for every number of the form 2, 4, p*^ or 2p'' with p 
an odd prime. On the other hand, Euler's theorem states that every invertible 
a £ Z/pZ* satisfies a'^^"^ = l[n]. Thus, for composite numbers n not possessing 
primitive roots, (p{n) is not a possible order of an invertible. We therefore 
use A(to), Carmichael's lambda function, the maximal order of an invertible 
element in the multiplicative group (Z/pZ*, x). See e.g. [161 [10113], for more 
details. Of course, A and (p coincide for 2, 4, p'' and 2p'^, for p and odd prime. 



^sw ox . com/ gmpl [maplesof t ■ com| [pari ■ math . u-bordeaux ■ f r [ [gap-system . org] 
[magma ■ maths ■ usyd . edu . au | 
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Then A(2'=) = 2^-^ Jqj. g > 3 -^^^^ ^j. ^he other cases, since j = 

Y[iPi ~ for distinct primes pi, we obtain this similar formula for A: 

A ^n^i^') ~ lcm{\{p'^^)}. Eventually, we also obtain this corollary of Euler's 
theorem: 

Corollary 4 Every invertible a within Z/pZ* satisfies a'*''") = l[n]. 

Proof. n = Y[Pi' for distinct primes pi. Then (p{p^^) divides X{n). This, 
together with Euler's theorem shows that a'^'-"-' = The Chinese theorem 

thus implies that the latter is also true modulo the product of the p^^ . □ 
This corollary shows that the order of any invertible must divide X{n). For n 

prime, the number of invertibles having order d|n — 1 is exactly (p{d) so that 
J2d\k = ^ for ^iT- ~ 1- We have the following analogue for n a composite 
number: 

Proposition 5 The number of invertibles having order d\X{n) is X^s^ 117=1 V'i'^j) 
for n — p\^ . . .p'^ and Sd = {(^i, ■ • ■ , d^) s.t. dj\Lp{p'^') and lcm{dj\ = d\. 

Proof. By the Chinese theorem, an element has order d if and only if the 1cm 
of its orders modulo the p'^' is d. Then there are exactly f{dj) elements of order 
dj modulo . □ 

Let us have a look of this behavior on an example: let n = 45 so that 
(p(45) = 6 X 4 = 24 and A(45) = 12. We thus know that any order modulo 9 
divides Lp{9) — 6 and that any order modulo 5 divides iy3(5) = 4. This gives 
the different orders of the 24 invertibles shown on table [1] It would be highly 
desirable to have tight bounds on those number of elements of a given order. 
Moreover, these bounds should be easily computable (e.g. not requiring some 
factorization !). In [5l[19], the following is proposed: 

Proposition 6 |3 Corollary 6.8] For n odd, the number of elements of order 
A(n) (primitive X~roots) is larger than ip{ip{n)). 

Now, this last result shows that actually quite a lot of elements are of maximal 
order modulo n. Using this fact, a modification of algorithm 1 can then produce 
with high probability an element of maximal order even though n is composite. 

6 Applications 

Of course, our generation can be applied to any application requiring the use of 
primitive roots. In this section we show the speed of our method compared to 
generation of primes with known factorization and propose a generalization of 
Miller- Rabin probabilistic primality test and of Davenport's strengthenings [7]. 
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Table 1: Elements of a given order modulo 45 



6.1 Faster pseudo random generators construction or key 
exchange 

The use of a generator and a big prime is the core of many cryptographic 
protocols. Among them are Blum-Micali pseudo-random generators W , Diffie- 
Hellman key exchange 0, etc. 

In this section we just compare the generation of primes with known factoriza- 
tion [1], so that primitive roots of primes with any given size are computable. 
The idea in [4] is to iteratively and randomly build primes so that the factoriza- 
tions of Pi — 1 are known. For cryptanalysis reasons their original method selects 
the primes and primitive roots bit by bit and is therefore quite slow. On figure 
|3]we then present also a third way, which is to generate the prime with known 
factorization as in [IJ, but then to generate the primitive root deterministically 
with our algorithm (since the factorization of p — 1 is known) . We compare this 
method with the following full-probabilistic way: 

1. By trial and error generate a probable prime (e.g. a prime passing several 
Miller- Rabin tests [18]). 

2. Generate a probable primitive root by Heuristic 2. 

We see on figure |4] that our method is faster and allows for the use of bigger 
primes/generators. 
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Prime size 

Figure 4: Blum-Micali primes with known factorization vs Industrial-strength 
primitive roots 

6.2 Probabilistic Lucas primality test 

The deterministic primaUty test of Lucas is actually the existence of primitive 
roots: 

Theorem 7 (Lucas) Let p > Q. If one can find an a > such that a^^^ = 
1 modp and a i ^1 modp, as soon as q divides p — 1, then p is prime. 

We propose here as a probabilistic primality test to try to build a primitive root. 
If one succeeds then the number is prime with high probability else it is either 
proven composite or composite with a high probability. 

Now for the complexity, we do not pretend to challenge Miller-Rabin test for 
speed ! Well, one often needs to perform several Miller-Rabin tests with distinct 
witnesses, so that the probability of being prime increases. Our idea is the 
following: since one tests several witnesses, why not use them as factors of our 
probable primitive root ! This idea can then be viewed as a generalization of 
Miller-Rabin: we not only test for orders of the form ^^^^ but also for each order 
of the form where g is a small prime factor of n — 1 . The effective complexity 
(save maybe from the partial factorization) will not suffer and the probability 
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can jump as soon as an element with very high order is generated. The algorithm 
is then a slight modification of algorithm [1] where we let F{B,Q) = 1 — (1 + 



Algorithm 2: Probabilistic Lucas primality test 
Input: n > 3, odd. 

Input: A failure probability < e < 1. 

Output: Whether n is prime and a certificate of primality, 
Output: or n is composite and a factor (or just a Fermat witness), 
Output: or n is prime with probability of error less than e, 
Output: or n is composite with probability of error less than e. 
begin 

Set P = 1, a = 1, = n - 1 and q ^ 2. 

2 

while Q > n's do 

Randomly choose a mod n. 

if gcd{a,n) 7^ 1 or gcd{a~^ — l,n) ^ {l;ri} or a"-^^ ^ l[n] or 
(q == 2 and n is not a strong pseudoprime to the base a) then 
|_ return n is composite. 

else if a 9 =1 mod n then 
Set P = P/q. 
if P <e then 

|_ return n is probably composite with error less than P. 

else 

- Set e to the greatest power of q dividing Q. 

- Set Q = 

- Set a — a X a 1" . 

- Set k = kU{q''}. 

- Refine B such that F{B, Q) 4e. 

- Find a new prime factor q oi Q with q < B, otherwise set 
L lq = Q- 

if Every q was prime then 

I return n is prime and (a, k) is a certificate. 
else 

|_ return n is probably prime with error less than F{B, q). 

end 



Remark 8 The exponentiations by can in practice be factorized in a "Lucas- 
tree" IMW- 

Remark 9 Algorithm [H is correct for the primes and most of the composite 
numbers. 

Proof. Correctness for prime numbers is the correctness of the pseudo primi- 
tive root generation. 
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Now for composite numbers: the idea is that first of all, only Carmichael num- 
bers will be able to pass the pseudo prime test several times. 
The 4e then follows since at least one a passed the strong pseudoprime test. 
This reduces the possible Carmichael numbers able to pass our test. Then, for 
most of the Carmichael numbers, A(n) divides n — 1 but, moreover, A(n) also 
divides ^^^^ for some factor of n — 1. Therefore, a i will always be one. If 
n is prime on the contrary, only ^ elements will have order a multiple of q. 

Now for the in the loop. The argument is the same as for the Pocklington 
theorem [6, Theorem 4.1.4] and the Brillhart, Lehmer and Selfridge theorem f6l 
Theorem 4.1.5]: let ?i — 1 = kQ and let p be a prime factor of n. The algorithm 
has found an a verifying a"^^ = 1 mod n. Hence, the order of a'^ mod p is a 
divisor of = k. Now, since gcd{a i — l,n) — 1 for each prime q dividing 
k, this order is not a proper divisor of k, so is equal to k. Hence, k must be a 
divisor of p — I ~ v{p)- We conclude that each prime factor of n must exceed 
k. From this, Pocklington's theorem states that if k is greater than ^/n^ n is 
prime. And then, Brillhart-Lehmer-Selfridgc theorem states that if k is in be- 
tween n's and then n must be prime or composite with exactly two prime 
factors [6^, Theorem 4.1.5]. But n has escaped our previous tests only if n is 
a Carmichael number. Fortunately, Carmichael numbers must have at least 3 
factors [T71 Proposition V.1.3]. Now, whenever Q is below na, k exceeds 
and then n must be prime otherwise n would have more than 3 factors each of 
those being greater than ns. □ 

Here is an example of Carmichael number, 1729. 1728 — 2^3"^, where 
A(1729) = 2^3^. Then ^ is either 864 or 576 both of which are divisible 
by 36 = A(1729). Therefore, our test will detect 1729 to be probably compos- 
ite with any probability of correctness. Figure [5] shows that this algorithm is 
highly competitive with repeated applications of GMP's strong pseudo prime 
test (i.e. with the same estimated probability of correctness). Depending on 
the success of the partial factorization, our test can even be faster (timing, on a 
PIV 2.4GHz, presented on figure [5] are the mean time between 4 distinct runs). 

Haplessly, some Carmichael numbers will still pass our test. The following 
results, sharpening |1H lemma 1], explains why: 

Theorem 10 Let n — ...p'^. Let q be a prime divisor of (p(n), and 
(/i, . . . , f^) be the maximal values for which q-^' divides ip{p1^). There are 

invertible elements of order divisible by q (i.e. for which a i ^1 mod n). 

Proof. By the Chinese remainder theorem, one can consider the moduli by 
p^^ separately. Suppose, without loss of generality, that pl^ is such that /i > 0. 
Otherwise all the fi are and the theorem is still correct. Consider a generator 
g of the invertibles modulo Pi^ . An element has q in its order if and only if its 
index with respect to g contains q^^ . There are exactly 1 such elements 
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Figure 5: Probabilistic Lucas vs GMP's Miller- Rabin for primes with probability 

< lo-*^ 

among the elements of Z/p^^Z. By the Chinese theorem, among the elements 
having their order divisible by q modulo n, we have then identified (p{n){l — ^) 
of them: the ones having their order modulo divisible by q. Now the others 
are among the <(2(?t-)(^) that remains. Just now consider those modulo P2^. 
If /2 == then we have not found any new element. Otherwise, 1 of 
them are of order divisible by q. Well, actually, in both cases, we can state that 
1 — of them are of order divisible by q. We have thus found some other 

elements: (p{n){-^){l — ^). This added to the previously found elements 
makes (p{n){l — ^jj^tj). Doing such a counting for each of the remaining p^' 
gives the announced formula. □ 
For instance, take a Carmichael number still passing our test whenever 
B < 1450: 37690903213 = 229 x 2243 x 73379. WeU, 37690903212 = 19 x 2^ x 
3x59xl451xl931and A(37690903213) = 19x22x3x59x1931. Then, Q will be 

1451 X 1931 and our algorithm will be able to find elements for which ^ 1 

mod n: those of which order is divisible by 1931. Unfortunately, there are 
quite a lot of them: (p(n)i|§ = 37489647840 w (1 - .00533962722683134975)n. 
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Thus, there are more than 5 chances over a thousand to choose an element a 
for which a "51x1931 ^ \ mod n. Even though this is much higher than ^ (if 
n was prime) , this probabihty will not be detected abnormal by our algorithm. 
Now, even if p — 1 is seldom smooth for p prime [21 , one can wonder if this is 
still the case for this special kind of Carmichael numbers . . . 

7 Conclusion 

We provide here a new very fast and efficient algorithm generating primitive 
roots. On the one hand, the algorithm has a polynomial time bit complexity 
when all existing algorithms where exponential. This is for instance illustrated 
when comparing it to existing software on figure O On the other hand, our 
algorithm is probabilistic in the sense that the answer might not be a primitive 
root. We have seen in this paper however, that the chances that an incorrect 
answer is given are less important than say "hardware malfunctions" . For this 
reason, we call our answers "Industrial-strength" primitive roots. 

Then, we propose a new probabilistic primality test using this primitive root 
generation. This test can be viewed as a generalization of Miller-Rabin's test 
to other small prime factors dividing n — 1 The test is then quantifying the 
information gained by finding elements of large order modulo n. When a given 
probability of correctness is desirable for the test, our algorithm is heuristically 
competitive with repeated applications of Miller-Rabin's. 
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